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ABSTRACT 


A hypothesis is made that the Galerkin Finite Element 
Method (GFEM) offers a viable option to the traditional 
Finite Difference Method (FDM) for numerical weather predic- 
tion. The shallow water Darotropic primitive equations are 
the forecast equations for all experiments. The hypothesis 
is tested by observing simple, analytic, atmospheric wave 
propagation on uniform and variable mesh grids. Second, a 
strongly forced solution simulating small scale nonlinear 
interactions is evaluated for both the GFEM and FDM. 
Finally, a variable, moving grid for a GFEM model is 
compared to a uniform, higher resolution GFEM model for a 
strong vortex ina mean flow. The G*EM shows a better 
propagation for simple atmospheric waves and better predic- 
tion to a forced nonlinear solution than the FDM model. A 
moving variable grid follows an area of strong gradients 


while not generating noise in the transition zone. 
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I. INTRODUCTION 


Proper simulation of atmospheric flow is the main 
objective of numerical weather prediction. The foundation 
of numerical weather prediction is a set of equations 
including the momentum, thermodynamic, moisture and conti- 
nuity equations for modeling the atmosphere. Through 
computer simulation, numerical weather prediction predicts a 
future state based on initial conditions describing the 
atmosphere. A successful forecast model must include small 
scale processes. Transports, conversions and exchanges of 
mass, momentum and energy occurring on the small scale 
represent important features which must be properly simula- 
ted. Feedbacks from the proper representation of these 
small features can sometimes markedly influence the larger 
scale solutions. Representation of the effects of small 
scale processes can be accomplished directly through 
increased spatial resolution. Continued increases in the 
spatial resolution will eventually allow the desired process 
to be resolved properly. However, a doubling of resolution 
generally requires an eight fold increase in computational 
effort. The value of increased spatial resolution must 
ultimately be measured by its contribution to the overall 
forecast. A level of confidence in the forecast must be 


achieved that the process is resolved near that particular 
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grid resolution. The grid resolution could be uniformly 
fine in the forecast domain. The computationa! effort for 
uniform fine mesh models could be far in excess of available 
resources. An alternative to a uniform fine mesh is a 
variable mesh where the fine mesh covers only regions of 
interest or high activity. 

The conventional numerical weather prediction forecast 
scheme, the finite difference method, approximates the 
partial differential equations with a truncated Taylor 
series. It has performed admirably when forecasting the 
larger scales of motion. Technological improvements in 
computing power coupled with better understanding of the 
atmosphere now allow the smaller scales to be fsrecist. 

An alternate method for numerical weather prediction, 
the GFEM MIO UH ES Ка partial differential 2quations 
while minimizing the error between the actual equations and 
their approximation. This best fit logicaliy leads one to 
the expectation that the GFEM will better model the smaller 
scales than the finite difference schemes. This research 
will demonstrate practical aspects of the GFEM tneoretically 
possible. 

In this dissertation, the GFEM will be evaluated to 
determine its potential to model atmospheric flow. Equiva- 
lent GFEM and FDM models will be utilized to compare the two 
metnods. Simulations of small scale processes explicitly 


resolved on uniform and variable grids will be investigated. 
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Rigorous experiments will demonstrate both GFEM and FDM 
responses for forcing near the grid length scale. Finally, 
a demonstration will be made of a GFEM moving variable mesh 
which moves with the small scale process or feature, and 
thereby allows the fine mesh to resolve the highly active 


region initially and throughout the forecast period. 
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ШЕ HYPOTHESIS 


Numerical weather prediction has steadily improved in 
the simulation of atmospheric flow. Large scale flows have 
been adequately represented by finite difference models for 
several decades. Galerkin-type formulations (Cullen, 1974b; 
Hinsman, 1975; Staniforth and Mitchell, 1978; Cullen and 
Hall, 1979; Staniforth and Daley, 1979; MacPherson and 
Aksel, 1980; and Sasaki and Reddy, 1980) have been shown to 
be competitive with finite difference models, but have not 
shown a marked improvement. Comparing the current opera- 
tional models at selected large computer centers, one finds 
that two are Galerkin and three are finite difference 
(Galerkin: National Meteorological Center and Canadian 
Meteorological Center; and finite difference: Fleet 
Numerical Oceanography Center, Air Force Global Weather 
Center and European Center for Medium Range Weather 
Forecasting. However, all centers have ongoing research 
with Galerkin models and there are indications that these 
will give better long range forecasts. The dichotomy arises 
because the Galerkin applications have not vindicated them- 
Selves with a marked increase in accuracy, but rather have 
Shown equivalent accuracies. Staniforth and Mitchel (1977) 
Stated that ultimately the best global/hemispheric models 


will be a spectral model. The spectral model is based on 
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the Galerkin procedure and the use of trigonometric func- 
tions is particularly appealing for hemispheric or global 
grids. However, where non-uniform grids are required, the 
GFEM is a more logical choice. 

Marked improvements in regional forecasts will not come 
from better large-scale models, but rather from models that 
can simulate smaller atmospheric features. The National 
Weather Service Limited Fine Mesh Model is an improvement 
over the hemispheric model, partly because it has higher 
resolution and can resolve smaller phenomena. These smaller 
features can also affect the large scale flow. Numerical 
meteorologists have realized this for years and have 
attempted to model these features. 

Mime nas the potential to increase efficiently the 
spatial resolution for the purpose of simulating accurately 
the small-scale processes. `f a variable mesh is to be 
employed, then it shouid be 2valuated by proper simulation 
of a small-scale feature. In previous research (Staniforth 
and Mitchell, 1978), the refined grid was tested by pro- 
pagating synoptic-scale waves into the finer grid and 
demonstrating that the wave could move into the finer grid 
without generation of significant noise. These conditions 
represent a prerequisite. If synoptic-scale waves cannot 
move freely into and through the variable grid, then advec- 


tion interactions will not occur properly in the fine grid. 
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However, one must also show that a smaller scale atmospheric 
feature can be properly resolved or developed in the fine 
mesh area. This will be a milestone. The efforts of 
Staniforth and Mitchell (1978) and Staniforth and Daley 
(1979) have stressed the movement of synoptic scale features 
into the fine mesh area but have not demonstrated improved 
resolution of smaller scale atmospheric features in the fine 
mesh area. 

The hypcthesis is that the GFEM isa viable option for 
numerical weather prediction when simulating atmospheric 
flow on variable grids. Three separate features of the GFEM 
will be explored. Each feature will establish the creden- 
tials of the GFEM as a viable option. Each feature is 
intimately related to the proper representation of a small- 
scale phenomenon. The cost effectiveness of the GFEM, when 
evaluated in these contexts, will provide a measure of the 


potential contribution of the method. 


A. FEATURES 
The following three specific features will be explored 
to support the hypothesis. 
las Variable Grid 
Investigations of a suitable alternative to the 
finite difference models for a variable grid will be 
performed. Two basic subdivisions are available--triangular 


or rectangular. Staniforth and Mitchell (1978) employed the 
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variable rectangular grid as shown in Figure 1. The grid 
contained some areas of finer resolution which were not in 
the verification area. Two points associated with having 
variable resolution in peripheral regions arise: 


(a) Unnecessary computational overhead is required; 
and 


(b) Undesirable phase changes occur as the wave 
propagates in the peripheral regions. 


Older (1981) and Woodward (1981) developed г transformation 
procedure to vary smoothly the resolution for a channel 
domain from a coarse toa fine area in a triangular subdivi- 
sion. A uniform equilateral subdivision is shown in Figure 
2. The use of triangles allows the increased resolution to 
occur only where desired and not in peripheral regions. 
Two possible choices of triangles include the right and the 
equilateral. Hinsman (1975) utilized equilateral triangles 
while Cullen (1974b) utilized near-equilateral triangles. Both 
reported excellent wave propagation. Kelley and Williams 
(1976) utilized right triangles and experienced very noisy 
solutions. Woodward (1981) duplicated Kelley and Williams 
effort with equilateral triangles and found a major 
reduction in the noise. 

The differing geometries and results thus far reported 
mandate a further review of triangular subdivisions. It is 
not obvious which subdivision is most suitable. A distinct 


advantage of the rectangular subdivision is that it allows 
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algorithms to be developed which take full advantage of 
vector processors. However, both subdivisions afford the 
luxury of obtaining a variable grid. Therefore, similar 
experiments performed on each type subdivision will point 
Out the advantages and disadvantages of each. 

2% Numerical Simulation of Physical Processes 

The variable grid chosen as the most suitable 

alternative will be tested while simulating a physical 
process to show that a small-scale atmospheric feature can 
be properly portrayed in the fine mesh area. As stated in 
Chapter II, the ability of variable grids to resolve small 
scale atmospheric features has not been rigorously tested. 
The ability to move synoptic scale features into a finer 
mesh is a prerequisite and must be shown. The main crux of 
the problem is what is happening in the fine mesh area. A 
proper scheme must resolve a small scale feature near the 
smallest grid length. Schoenstadt (1980) and Williams 
(1981) indicated that the most responsive schemes were either 
a staggered finite element grid with primitive variables or 
an unstaggered finite element grid with vorticity/divergence 
formulation. This research will utilize the latter. 

The simulated process will be that of a mass source 
analogous to that found in the upper atmosphere above a 
hurricane or on the leading side of a strong trough. The 


source will appear in the continuity equation. A known wave 
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form will be assumed and the required source to support the 
wave form will be analytically derived. The expression for 
the wave form will allow control of the scale of the process 
so that a near-grid scale phenomena can be simulated. 
За Moving Grid 

The ability to move the fine grid so that it 
remains centered on an atmospheric circulation will also be 
demonstrated. This capability will further enhance the 
applicability of the GFEM for atmospheric simulation. If it 
is shown that small-scale forcings are proper!y reflected in 
the flow, and the grid can be moved with the forcing 
disturbance, then the GFEM has great potential for 


atmospheric prediction. 


24 





INL LL IHPSBEVELOPMENT ОР THE СРЕМ 


The hypothesis in Chapter II proclaims the GFEM as a 
viable option for simulating atmospheric flow. An under- 
standing of the history and principles of the GFEM will give 
further insight for the expected improved performance as 
compared to the conventional FDM. 

The origin of the GFEM can be traced to the seventeenth 
century work by Leonhard Euler. His work established the 
branch of mathematics known as calculus of variation. He 
recognized traat there existed a partial differential equa- 
tion (PDE) associated with the minimization of a functional. 
This PDE has been appropriately named the Euler-Lagrange 
equation. Solving the PDE was equivalent to determining a 
stationary value of the functional. In the ес een 
century, Lord Rayleigh furthered the variational calculus by 
representing the dependent variable as a mathematical 
expression, typically as a power series. The stationarity 
of the solution allowed determination of the unknown coeffi- 
cients of the power series. His work was generalized in the 
early twentieth century by W. Ritz and the procedure has 
Since been referred to as the Rayleigh-Ritz method. A 
limitation in the solution by the Rayleigh-Ritz method is 
the fullness of the matrix to be inverted. Shortly after 


Ritz completed his work, Galerkin developed a procedure 
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using weighted residuals which had a wider application than 
classical variational calculus. The residual is made ortho- 
gonal to a function, called the test function, and thereby is 
minimized. While the Rayleigh-Ritz procedure applies to the 
functional of an associated Euler-Lagrange equation, the 
Galerkin procedure pertains to any PDE. The methods produce 
identical results when applied to equivalent extremum 
problems. 

The Galerkin procedure was unknowingly well establisned 
by the end of World War II. The rapid technological 
advances made during and after the war, coupled with tne 
emergence of the computer, led the aircraft industry to 
develop new numerical procedures for solving stress problems 
in aircraft design. A successful technique was developed 
and mathematicians realized, long after the fact, that the 
Galerkin procedure was being utilized. 

The GFEM is a commonly used subset of the set of 
Galerkin procedures. After subdivision into a set of ele- 
ments, the domain resembles a completed jigsaw puzzle, and 
hence leads to the terminology "finite elements." Figures 1 
and 2 are samples of such subdivisions. The dependent 
variables of the PDE are represented as a linear combination 
of known functions, which are usually low order polynomials. 
The same function is employed as the test function. When 
Substituted into the PDE, these approximations leave a 


residual. Minimizing the residual completes the procedure. 
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The functions are globally zero except where the dependent 
variables are defined near each individual node. The pro- 
duct of functions remains zero except where the representing 
and test functions are both non-zero. This makes the method 
attractive for computer implementation. It removes the 
matrix "fullness" problem found with the Rayleigh-Ritz 
approach. The minimization process lies at the heart of the 
expected improved performance. Each term of the PDE has 
been simultanecusly approximated and the error in those 


approximations has been minimized. 
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IV. MODEL DESCRIPTION 


Two different barotropic shallow-water models will be 
employed to test the hypothesis in Chapter II. Each model 
will have the identical domain, boundary and initial 
conditions. The difference will lay with the partial 
differential equation approximation. In one model, a 
Galerkin approximation to the partial differential equation 
will be used, while the other model will use a finite 
difference approximation. Two versions of the Galerkin 
model will be evaluated. The first version will have trian- 
gular subdivisions and basis functions, while the other will 


have rectangular subdivisions and basis functions. 


A. GALERKIN FINITE ELEMENT MODEL 
The system of equations referred to as the shallow- 
water equations consists of three equations with three 


forecast variables 9, u and v. The equations are written 


3$ , 2d 4 yb 4 (QU OS 4-1 
и озуу 0 a? 
au , 3U, 3U аф 4-2 
3t ^ "ax * V3 et 0 ве 
+ еШ + EM + fu + 39 . 0 (4-3) 
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Here $ 1$ the geopotential height, 
u is the east/west component of the wind, 
v is the north/south component of the wind, and 
f is the Coriolis parameter 
By expanding » into a mean (6) and a deviation (4') the 


equations can be written 


2%" E ыд. шы, in 

хэ Ээ даан) 3-1 у 0997.) 80 ШЕ? 
ди + 9%! + ak - vQ = 0 (4-5) 
9% ð X 9X 

ӘУ yk : 

ло) “зу? uQ = 0 (4-6) 


here Dis the divergence 

K is the kinetic energy (per unit mass), and 

Q is the absolute vorticity. 
The primes will be dropped for the rest of the paper for 
сщагтъу. 

Cullen and Hall (1979) showed that the accuracy of the 
GFEM solution was better for the vorticity-divergence 
formulation of the shallow-water equations than for an 
increase in resolution with the primitive formulation. 
Williams and Schoenstadt (1980) noted that staggered 
variable formulation of the primitive equations and the 
unstaggered vorticity-divergence formulation gave the best 
treatment of geostrophic adjustment for small-scale 


features. 
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The vorticity-divergence form of the shallow-water 
equations also allows the use of a semi-implicit time 
scheme. This scheme artificially slows the propagation 
Speed of the fastest gravity waves, which allows a much 
larger time step than one could expect for a normal Courant- 
Fredrich-Lewy (CFL) stability criterion. This scheme thus 
offsets some of the extra computational expense required to 
solve the system of equations assembled at each time step. 


The vorticity/divergence form of the equations is 





9 9 9 

22 * 9D + 27049) Е zy Vo) 2 (4-7) 

BC os 9, = - 

E. 00 (vq) «0 (4-8) 

90 2 2 д 9 " 

bu "y v^K e zx Md) + sy (ua) 2018 (4-9) 
Here vé is the Laplacian operator апа с is the С = A 


relative vorticity. 
The velocity can be writzen a5 the sum of the rotational and 


irrotational components as 


Инеге у = KiVy and I, = УХ 


Y 


The equations can be rewritten using 


u ду 
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D = 3x 7 3 X 
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resulting in 


3 2 3 3 
22. Emm ay 99) DS (ve) (4-10) 
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Equation (4-12) can be further manipulated 


75128 99) 4 24040) - 2E - 320008) + 2) (4-13) 


The domain of integration is a channel with east-west 
periodicity. The boundary condition at a wall is 
Yo N= 0 
where N is an 1) pointing normal vector. Along the 


northern and southern walls, the v component is equal to 


zero, so that the v equation of motion (4-3) reduces to 
шю _ 
Зу Ти 

The zonal and meridional components of the wind can be 


written as non-divergent and irrotational components as 


цаг REMO. SX 
ду ð X 


and 


121201 
aX ду 


Then, along the north/south walls where v equals zero, the 


boundary condition is 


(4-14) 
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The above condition is imposed by setting 

p = a constant (4-15) 
when solving the vorticity equation and 

= = 0 (4-16) 
when solving the divergence equation. This is an overspe- 
cification but (4-14) would be difficult to apply. The 
initial conditions which will be presented later are 
specifically selected to satisfy (4-15) and (4-15). 

The semi-implicit scheme is implemented by evaluating 
all the terms on the left hand side of the equations as an 
average at time levels (t + At) and (t - at) or with a 
centered time difference as appropriate. All the terms on 
the right hand side are evaluated at time level t. The 


equations become 


2 + o(0°y(ttat) + 1 = - =5-(00) - 371У2) (4-17) 

7 A . ` = (uQ) т зид) (4-18) 

цаг: + o(ttat) + H(t-At)) = s; (v) : = а 52.0000) Т 35) 
(4-19) 


The divergence equation (4-19) can be solved for x at 


(t * ^t)and substituted into the equation (4-17) to yield 
" E 
A Oe $59 - ша), 55) 


4-20) 
sl), 1 2 5 a ( 
(АДЕ ЫС) Kitu) + syltv)) 
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where the overbar denotes an average of t + At and t - ^t. 
The system of three equations in three unknowns has been 
reduced to two Poisson equations and one Helmholtz equation 
to be solved at each time step. 

The solution procedure involves solving the » equation 
(4-20) for a new %. The divergence equation (4-19) is then 
solved for Ф + 2% and by subtraction for 2%. Finally, the 
vorticity equation (4-18) is solved for =. There are two 
options as to whicn variables will be history-carrying: 
either ф, № апа х ог ġġ, u and v. The choice was made to use 


$, u and v as history-carrying variables. They are updated 


after each time step following Staniforth and Mitchell 


977) Бу 
e(ttat) = 26 - ó(t-At) (4-21) 
u(ttat) = 2at (= = ° m "x + u(t-at) (4-22) 
v( test) = 24% (5 E n E "d + v(t-At) (4-23) 


Implementation of the GFEM is accomplished as described 
in Chapter III. The north and south latitudes are input 
parameters, and the north/south distance is subdivided into 
N equal parts, where N is also an input parameter. The 
east/west increment is a function of the north/south incre- 


ment, as explained in the section describing the individual 
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experiments. The Coriolis parameter will be a constant 
equal to the mid-channel value. 

The triangular and rectangular domain subdivisions will 
be demonstrated in the experiment sections. An appropriate 
approximating function for the rectangular subdivision is a 
bilinear function. In either case, the forecast equations 


in Galerkin form are 


QAM 
м 9 > ША 
[0555 - u O Y а қр 
ф. 1 
A үр _ 
v ДАРЬ Л xtt- At)V, V. (4-24) 
2 9) _д_ _9_ 
| = В Nr (uQ); ХН (v0) ; o T 
2 9. " = д 9K 
aK (4-26) 
- бу об) * y V0 


Here the integral sign implies an area integral over the 
domain, the j subscript denotes Einstein summation for the 
dependent variables and the i subscript is the ith nodal 


equation. 
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The strength of the GFEM is that the spatial deriva- 
tives of the dependent variables become derivatives of the 
known approximating function and, therefore, are known 
exactly at each node. When using piecewise continuous 
linear functions, the first derivatives are piecewise 
discontinuous and second derivatives are not defined. 
Therefore, the second derivatives must be handled in a 
different fashion, 


% 22 92 
f" V оу, —3( Vi) 


ду 


3 D» fav. av. 
EE EL TUE ) ё эх! 
av. av. 
bx Sx ji | әрі ax | 


9 s 202 ӘУ. V ду au. 
СЕЕ у. әу 
/ 2411 ayvay) |) 7 Ei | 


= av. ji ^ MEC 
$ è ay" : 77 зу! 


where the $ implies a iine integral alóng the east/west 
or north/south boundaries. The east/west line integrals 
are zero because of periodicity,but the north/south line 
integrals add additional terms to the equations and are 
satisfied by the boundary conditions. The Laplacian terms 


become 
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This integration-by-parts procedure allows the use of a 
linear function while approximating second derivatives. The 
final form of the forecast equationsis 

ЖАДА 
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(4-27) 
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Here the Helmholtz equation (4-27) has utilized 


pr г. JU , 9V 

v^x * DIVERGENCE - 3T * зу 
and 

29 -. F 

3y ur, 


along the walls. 

The line integral along the north and south boundary 
has been dropped from the vorticity ecuation (4-28), since 
the value of y on the boundaries is a constant and therefore 
the value of 25. is zero. The line integral along the north 
and south boundary has been dropped from the divergence 
equation (4-29), because the value of the normal derivative 
along the north/south walls is zero, es stipulated by the 
boundary conditions. The initial conditions will also 
satisfy the condition that the normal derivative of x is 


zero along the north/south walls. 


Br PERITE DIFFERENCE MODEL 

The comparison model to the GFEM model is an adaptation 
of the staggered, primitive-equation model as described in 
Section 7-4 of Haltiner and Williams (1980). The equations 
are in the flux form and employ Scheme C as shown in Section 
7-3 of Haltiner and Williams (1980). The baroclinic model 
described there has been coded to perform in a barotropic 


mode. 
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Since the GFEM is written in differentiated form, 
vorticity and divergence have been diagnosed from the primi- 
tive variables in the finite difference model to allow 


appropriate comparisons. 


е. NUMERICAL METHODS 

Although no attempt is made to optimize computational 
efficiency in this research model, the GFEM must utilize 
efficient numerical techniques to be considered a viable 
option for numerical weather prediction. Various solution 
procedures are available for the GFEM system of equations. 
Successive over-relaxation was employed during the early 
research stages. Several problems arose which indicated a 
need for a better solution procedure. The successive over- 
relaxation procedure resulted in a bias if the sweeping 
during each pass was in the same direction. Alternating the 
direction alleviated that particular problem, but the number 
of passes required to achieve a desired level of accuracy 
became prohibitive. A second approach employed a direct 
SOlver using a Gaussian elimination procedure. The matrices 
from the Galerkin procedure are decomposed into upper and 
lower block tri-diagonal matrices. A preprocessing, repre- 
senting the forward substitution stage, can be done once. 
The back substitution must be done each time a solution is 
desired. The coefficients necessary for this back substitu- 


tion are stored in an efficient manner. The particular 
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algorithm is called a “skyline” solver and is described in 
detail in Bathe and Wilson (1976). Skyline refers to the 
compact method of storing only those coefficients required. 
This method has the desired level of accuracy and a high 
degree of computational efficiency. Staniforthn and Mitchell 
(1977) employed another direct solver method called, the 
conjugate- gradient method. This method is very computa- 
tionally efficient, even when the non-constant coefficient 
Helmholtz equation is involved. The theoretical operation 
count indicates that for very large domains with many’ 
degrees of freedom, the conjugate-gradient method wi 1 be 
much better than the other direct solver mentioned adove. 
The Gauss elimination procedure was utilized instead of the 
conjugate-gradient method for reasons of expediency. 
Numerical integration of the three forecast equations 
involves solving first a Helmholtz equation for о, second a 
Poisson problem for y and finally a Poisson problem for x. 
The boundary conditions and the equations are well posed for 
the first two equations. However, the normal derivative of 
x is equal to zero for the last Poisson equation and 
requires special attention. The solution plus a constant is 
also a solution. Since only derivatives of y are required, 
the shape of the field is much more important than the 
actual value of y. To avoid computer round off errors, the 
average value of the terms on the right hand side of (4-29) 


1s removed at each time step. This sets the constant to 
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zero and does not allow the solution to grow by a constant 
each iteration. 

The current tecnnology of large mainframe computers 
incorporates vector nrocessing. Utilization of this 
ability can be acnieved when the system of equations can be 
уесіогігей. Тһе СГЕМ 15 derived in vector formulation and 


is easily vectorized. 


D. INITIAL CONDITIONS 
l. simple Atmospheric Waves 
The initia! conditions must allow a relative 
amount of control for the desired input parameters as well 
as satisfy the beundary conditions. The forecast model 
history-carrying variables are d, u and v. The analytic 


expression for the streamfunction, у, is 


| 552078 22 дтһХ г 
U = A sin E sin eS = Uly-y га) + = (4-30) 


where A = ampiitude of the perturbation, 
W= width of the channel, 
L = length of the channel, 
n = wave number, 
U = mean flow speed, and 


y ume Wddle post of the channel 
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The first term in the expression is the perturbation. 
The second term represents the north/south slope necessary 
to support amean flow of U. The third term is the mean 
depth term. The geopotential height, o, is related geostro- 


phically to the streamfunction, y, by 


By use of a trigonometric identity, the streamfunction is 


written as 


y = 50 = cos ZZ) s | #rnx - Uly-y Hs 


$ 
mid. т, 


The vorticity is given by 
pe- Фа ГА 51Па2Х COS 2a4y - a2? A sinagx 


2 (4-33) 
+ а? 5 Sinagx совидату 


where 
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The divergence is determined through a linearized form of 
the quasi-geostrophic divergence equation from Chapter 3-2 


of Haltiner and Williams (1980), in the form 


a 111. (4-34) 
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The right-hand side of this equation becomes 


¡e 


4237] — (4-35) 


ГЬ БЕ. 
COS а2Х Г(-202а1“Аа2 +004227) соѕ дауу зайн 


Assuming a solution form for divergence on the left-hand 


side 
DISCOS aox(C41cos Фау + Co) 


{һе constants С, апа Co become 
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То - А Ч 
Ср = - (2 Aag + -p Шад ?2)/ (да + во + 6) 
апа 
2 
f 3 2 f 
Ё G A _0 
Co = ple, 7/lay +32) 


Using 72, = D md assuming a solution form for x as 

x = cosaox(C3cos2a4y + Сд) (4-36) 
then C4. and C4, become 

Cz = - C1/(4a4* * a2?) 
and 


C4 * - Cg/ag* 
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Observe that y equals a constant and № equals zero along 
the north/south walls. With the expressions for stream- 
function, (4-32), anc velocity potential, (4-36), appropriate 


formulae for u and v can be derived. 


u = U - sin aax(Acı sin2aıy + ao2(C4 cos Фау + Ca) (4-37) 


and 


A 
у = с05 а2Х(825(1-С08 2а1у) - 2C3a,sin 2a,y) (4-38) 


and 6 is given by equations (4-31) and (4-32). 
2. Source Term 

The addition of a source term into the continuity 
equation wiii test the resolvability of the GFEM model for a 
source. The source is constructed to represent a small- 
scale meteorological phenomenon to provide a measure of the 
response of the GFEM model to forcing near the smallest grid 
scale. The source initial conditions must be able to 
describe a smail-scale feature embedded in a large synoptic 
Situation. The initial conditions are derived by choosing a 
desired solution and then back substituting to derive the 
form of the source required. This is a unique approach when 
adding a source term. First, a source, S, is added to the 


continuity equation in the form of 


Re Puce er Se 


at ЭХ 3 y Sa (4-39) 
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Following the scale analysis of Chapter 3-2 in Haltiner and 
Williams (1980),a new quasi-geostrophic potential vorticity 


equation is formed with the source term included, 


2 2 f З 
Ы v syl - “B) + FS= 0 


Solving for S leaves 


2 2 
Шоо _ әф 9 221 123 9 
S = (st y ох ra Eu Ë #55 š a) ша 
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The technique for testing the source term is to pick a 
streamfunction and then solve for the source term which will 
Satisfy that streamfunction. The particular streamfunction 
evaluated will be a steady state solution. The source in 


this case is 


2 2 2 2 
= — ох 2 9-р, 39 9(9 9 ,9 Ú 
P a p Е Bu x T 7 +) (4-41) 
| ду ðX ду 
The particular form of the streamfunction is 
1 2 тх 
x 2 “отл ү) 
=Z a ы = 1 ту R L = 
y И (у йг) Asin Wes (4-42) 


where R is a constant. The use of the exponential allows a 
variation in the scale whereby a small scale feature can be 
embedded in a larger flow pattern. Substitution of this 


particular y into the expression for the source yields 
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E. NONLINEAR EFFECTS 

Cullen (1974a) discussed the requirement to simulate 
the cascade of energy to unresolved wavelengtns. To avoid 
accumulation of energy in the snortest resolvable wave- 
lengths, he developed a two-step method for computing the 
nonlinear terms. Implementation of the two-step method can 
be accomplished efficiently. In the formal development of 
the model equations, there are several nonlinear terms. The 
two-step method involves projecting both variables in a 


nonlinear term into Galerkin space and then performing the 
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multiplication. This procedure insures the best approxima- 
sion to that nonlinear term in a weighted-mean sense. 
Additionally, the computational effort for the two-step 
method is much less than the one-step method which involves 


computation of all the possible interaction coefficients. 


г. STABILITY CRITERION 


The stability criterion for the forecast equations is 


fhis is more restrictive than a normal Courant-Fredrich-Lewy 
CFL) condition. It is made up of the normal CFL criterion 
based on advection plus another effec: due to rotation. 
‘inertial motion). It provides an upper boundary to the 
inaximum allowable time step when using a semi-implicit 
scheme. In this research, time steps in excess of an hour 
are easily used. Care must be taken not to exceed the 
stability criterion when larger mean flows are experienced 
or when moving to higher latitudes (larger f). 

A more detailed discussion of the derivation of the 


Stability criterion can be found in Appendix A. 
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oe eee! SAND RESULTS 


Experiments will be performed involving the features in 
Chapter II to test the hypothesis. Four separate experi- 
ments will be performed. The first experiment wii! 
establish the overall performance characteristics of three 
specified models. The next three experiments wii! address 
each specific feature supporting the hypothesis. 

Standards of comparison must be established for each 
experiment. The basic comparison will be achieved through 
harmonic analysis of the initial and forecast fieids. The 
harmonic analysis is a double Fourier decomposition. The 
harmonic analysis requires that the data be on a uniform 
grid. Because of the triangular and variable grids, an 
interpolation is necessary. A fifth degree polynomial 
surface fitting routine from the International Mathematics 
and Statistics Library (IMSL) was employed to interpolate 
data on the non-uniform grids to a uniform spacing. First, 
a harmonic analysis is performed along columns of constant x 
to obtain the y structure of the initial conditions. Then a 
harmonic analysis of the amplitudes of each y wave number 
along rows gives a true double harmonic analysis. This 
double transform approach E ry due to the initial 
conditions. The use of sin pe in the y direction is the 


Same as using (1.0 - cos Au), The boundary conditions are 
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satisfied by a sine wave in the y direction. Therefore, an 
infinite sine series is needed to represent the cosine 
function. 

The theoretical! movement of each wave is obtained from 
an analysis of the shallow-water equations. The analysis in 
Section 2-6 of Haltiner and Williams (1980) predicts the 


wave speed as 


where 
C equals wave speed, 
H equals height, and 


u X direction wave number. 


By generalizing the streamfunction to a two-dimensional wave 


form 
p = A sin M cos u(x-ct) (5-2) 


and after furtner manipulation, the wave speed becomes 


C = D/U + FE + AD) 


where 
ф equals average geopotential height, 
ux equals x direction wave number, and 


у equals y direction wave number. 
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A. EXPERIMENT 1 

Establishing wave propagation characteristics on a 
uniform grid will illuminate the properties of each model. 
Three models will be compared: 

1% A triangular subdivision GFEM model; 

2% A rectangular subdivision GFEM model; and 

34 An equivalent finite difference model. 
The GFEM models will use the differentiated form of the 
shallow-water equations with a semi-implicit scheme. The 
finite difference model uses the undifferentiated form on a 
Staggered grid with a centered (leapfrog) time scheme. The 
triangular subdivision uses equilateral triangles. The 
relationship between base and height for an equilateral! 


triangle is 
BASE = 2 x HEIGHT/Y3 


The x distance therefore has the same relationship to the у 


distance 
AX = 2 ay/ "3 


To perform comparisons of similar models, the same x, y 
relationship is retained for the rectangular and finite 
difference models. 

The basic difference between the rectangular and 
triangular models will be in the approximating polynomials. 


The rectangular polynomials are bilinear while the 
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triangular polynomials are linear. The higher order 
polynomials shouid provide an increase in accuracy. During 
the integration process, many integrals require evaluation. 
Numerical quadrature could be utilized, however, a more 
efficient method is available through the use of natural 
coordinates. Quadrature with no error can be accomplished 
by formula with this method. A detailed description of the 
natural coordinate method is delineated in Appendix B for 
both triangles and rectangles. 

Diagrams for the grids for the triangular, rectangular 
and finite difference models are shown in Figures 2, 3 and 4 
respectively. ihe domain is 4896.0 km in the y direction 
and 5653.0 km in the x direction. Each model has 12 incre- 
ments in the x and y directions. This gives 156 degrees of 
freedom. All tests will be for a 48 hour time integration, 
a mean flow of 10 m/s and a mean depth of 1000 meters. A 
perturbation of 1 m/s will be added to the geopotential 
field, which includes the mean height plus the north/south 
slope required for the 10 m/s mean flow. This small pertur- 
bation will focus on the linear aspects of the model's 
capabilities. Experiments described below will evaluate 
larger perturbations, and hence the nonlinear interactions. 


Models are evaluated for wave numbers 1, 2, 3 and 4. 
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Rectangular uniform subdivision for a channel in 


Cartesian coordinates. 
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mig 4. Finite difference uniform grid for a channel in 
Cartesian coordinates. 
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Bc RESULTS 1 

The propagation characteristics of the GFEM models and 
the finite difference model on a uniform grid will be 
evaluated. Graphical displays will highlight synoptic 
Characteristics. Harmonic analysis will show phase 
propagation and wave amplitude distributions. 

The initial and 48-h forecasts for the conditions 
Stated in experiment 1 with wave number one on a triangular 
grid are shown in Figures 5 and 6. The southeast/northwest 
Orientation in the forecast fields is due to the sine wave- 
form in the y direction. This waveform with model boundary 
conditions requires an infinite number of y modes. 
Fortunately, the amplitudes of the higher wave numbers are 
small. Consequently, there are only a few important modes. 
The multiple modes cause a skewness in the solutions for all 
three models since the y-modes have different phase speeds. 
Additionally, an analytic expression for the solution 
requires an infinite sum of the modes. For this reason, no 
true solution has been computed to compare with the GFEM and 
FDM modeis. 

The divergence equation is the most sensitive equation 
of the three forecast equations. A mean depth of 1000 
meters was chosen because the divergence is large. There- 
fore, a critical evaluation of the sensitivity of the GFEM 


model to this important meteorological parameter can be 
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Initial conditions for the GFEM model with a 
triangular subdivision and y 75 Е опе. 
Contour intervals are 600 m Bor < potential 
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vVorticity and .2x10-/s-1 for 2210, Nodal 
points are denoted E an x. 
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performed. A divergence field which has propagated 
uniformly with the vorticity field is shown in Figure 6. 
Harmonic analysis (not shown) indicates divergence ampli- 
tudes which initially adjust and then oscillate slightly 
about a mean value of .1270x1077 s-l. The phases of all 
forecast fields are uniform. Synoptically, the forecast 
fields $ u, and v indicate a smoothly resolved atmospheric 
wave. 

A rectangular subdivision 48-h forecast for wave number 
one is shown in Figure 7. Comparisons between Figures 6 and 
7] show little if any difference which indicates that the 
triangular and rectangular subdivisions give comparable 
results. Harmonic analysis of the divergence fields 
indicates a small oscillation about a mean similar to that 
observed for the triangular subdivision. The 48-h forecast 
from the finite difference model for wave number one is 
shown in Figure 8. The equivalent finite difference model 
also has 12 increments in each direction. Comparisons 
between Figures 6,7 and 8 show that the finitedifference 
model vorticity field lags the GFEM models. Additionally, 
the divergence field has lost its areal extent and developed 
a strong gradient of divergence. 

In each of the wave number one tests there were l2 grid 
points representing the wave. Good propagation should be 


expected from all three of the models for a 12-increment 
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wave. As the wave number is increased with the same degrees 
of freedom, the wave resolution will be decreased, and 
clearer comparisons can be made. 

Figures 9 and 10 are the initial and 48-h GFEM model 
forecasts for wave number 2, triangular subdivision. 
Forecasts for the rectangular subdivision and the FDM model 
are shown in Figures 11 and 12. The triangular and 
rectangular subdivision forecasts are similar. However, the 
FDM model forecast displays additional divergence and 
vorticity centers near the boundaries. Figures 13 through 
16 area series similar to Figures 9 through 12 except for 
wave number three. These waves are represented by 4 grid 
points. There are indications in Figure 14 that energy is 
appearing in wave numbers other than wave number 3. 
Comparison of Figures 14 and 15 show that the triangular 
forecast contains more noise than the rectangular forecast. 
Figure 16 shows a divergence field from the FDM mode! which 
is very different from Figures 14 and 15. There is a 
north/south asymmetry in the u and divergence fields. The u 
cells in the northern regions appear to be expanding while 
the southern cells are decreasing. 

Table lis a composite of the harmonic analysis for 
the triangular, rectangular and finite difference wave 


number 3 case (v component amplitudes). The initial and 
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48-h v component amplitudes for wave number 3 in the x 
direction and the first six possible y modes are tabuiated 
below. 
TABLE 1 
Harmonic analysis (v component m/s) - wave number 3 
Triangles Rectangles Finite Differerce 


Wave Initial 48-Н Initial 48-H Initial 48-H 


1 2.12 2.14 2.20 2222 2-2) 2220 
2 .01 .04 .01 .02 .01 208 
3 .41 .43 .44 .44 .44 .42 
4 .00 .00 „00 „00 „00 „СО 
S 21015) .06 .06 .06 .06 EE 
6 .00 .00 .00 . 00 „00 > 


fhe minor differences in the initial state amplitudes 
are due to the differences between the triangular and 
rectangular approximation for the initial conditions. 
However, the 48-h amplitudes show that the rectangular 
version has slightly better preserved the initial state 
amplitudes than either the triangular or finite difference 
models. The finite difference model has transferred more 
amplitude to the other modes, especially in the second y 
mode. 

Figures 17 and 18 are the initial and 48-h GFEM model 
triangular subdivision forecasts for wave number 4. Wave 


number 4 is a severe test for this model. The energy 
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transfer to other modes is readily apparent. Additionally, 
both the u component of the wind and the divergence show 
considerable noise at the boundaries. Figure 19 is the 43-h 
forecast for the rectangular model. No energy transfer nor 
noise at the boundaries is apparent in Figure 19. Figure 20 
is the 48-h forecasts for the finite difference model. The 
asymmetry in the wave number three finite difference case 

is magnified for wave number four. The divergence field is 
heavily asymmetric. The u cells have developed a dipole 
from each original center. Comparisons of Figures 18, 19 
and 20 show the superiority of the rectangular forecast. 
While the triangular version is noisy, it still has 
maintained the main features better than the finite 
difference model. 


Table 2 is similar to Table 1 except for the wave number 


4 case. 
TABLE 2 
Harmonic analysis (v component m/s) - wave number 4 
Triangles Rectangles Finite Difference 


Wave Initial 48 -H Initial 48 -H Initial 48 -H 


1 2.68 2.60 2.94 2.96 2.94 Е „96 
2 .00 .03 .01 ‚10 ‚01 „05 
3 „51 „59 „58 159 EG .60 
4 .00 .00 400 .04 .00 US 
5 405 .08 .08 ‚08 .08 .08 
6 . 00 .00 .00 .00 .00 207 
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The harmonic analysis shows that the triangular version 
transfers more energy in the fifth y mode then does either 
other model. However, the rectangular model transfers more 
energy in the second y mode than does either other model. 


Table 3 compares the divergence for the wave number 4 case. 


TABLE 3 


Harmonic analysis (divergence 5-1) - wave number 4 
(All numbers are scaled 10 to the minus 10.) 


Triangles Rectangles Finite Difference 


Wave Initial 48 -H Initial 48-H Initial 48-H 


1 2682.0 1 "0" “292920 240122 290170 4462.0 
2 „0 140.8 5 AO 447.4 5934.0 
9 521.3 JT JS 440.7 545.5 2314.0 
4 .0 20.5 .0 81.2 „0 323.5 
5 62.6 22025 94.4 DERG 68.7 1382.0 
6 „0 10.5 0 8.4 21) 882.0 


The finite difference model has greatiy amplified the 
initial divergence. A choice between the triangular or 
rectangular grids version based only on these results would 
be purely subjective. However, the response of both grids 
at this high wave number in terms of divergence is highly 


superior to the FDM. 
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Figure 21 is the phase propagation diagram for the 
three models obtained from the double harmonic analysis. In 
al! cases, the percentage propagation is for the first y 
mode and the appropriate x mode. The phase propagation 
diagram indicates that the triangular and rectangular models 
are comparable for the low wave numbers and diverge slightly 
at higher wave numbers. In all cases the GFEM models have 
better wave propagation than the FDM model. 

The analysis performed on the uniform grids has shown 
that each GFEM model performs in a highly satisfactory 
fashion. There are minor differences in the phase 
propagation, while the rectangular version has an advantage 
in the control of energy transfers which manifest themselves 
in the divergence. The finite difference model performs as 
well as the GFEM models for the long waves. However, when 
forecasting the shorter waves, the finite difference model 


does not perform as well as either GFEM model. 


03 EXPERIMENT 2 

This experiment will investigate two options for a 
variable grid. Atmospheric wave propagation on variable 
triangular and rectangular grids will be the test vehicle. 
Since this FDM model has no capability on a variable grid, 
it will not be evaluated during this experiment. The 
comparisons will be between a rectangular and triangular 


GFEM model. 
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Phase propagation diagram for the triangular, 
rectangular and finite difference models from 
Experiment 1. 





There are many criteria for determination of a suitable 
variable grid. Both rectangles and triangles allow a 
successful implementation of a variable grid. One 
criterion is the accuracy of resolving the atmospheric waves 
that move through the variable grid. A second criterion is 
the increase in resolution. Another consideration is the 
ease with which the grid can be refined. 

Cullen (1979) stressed the importance of a smooth 
Variation in element size when increasing the resolution 
from coarse to fine. Older (1981) developed a technique for 
transforming a uniform triangular grid into a variable grid 
when working with a GFEM model. He demonstrated that a 
smoothiy varyirg grid allowed atmospheric waves to propagate 
through the transition zone with much less noise generation 
than an abruptly varying grid. Some finite difference 
models with nested grids, for example Harrison (1973), have 
abruptly varying resolution. A generalization of Older's 
techniques allows a similar transformation for a rectangular 
grid. 

The test cases for experiment 2 will be fora simple 
atmospheric wave. À mean depth of 5 km, a mean flow of 
10 m/s, wave number one and a perturbation of 1.0 m/s are 


employed. 
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D. RESMETS 2 
The ability to refine a grid should be easily achieved. 
The method of transforming from a uniform grid to smoothly 
varying should allow one to choose not only what degree of 
refinement but also where the refinement occurs. Older 


(1981) developed a transformation method based on 
X = X + A cos Kx 


where 
X = the transformed grid 
x = the original grid 
A = a constant 
к = 21/1. 


Woodward (1981) modified the transformation to include a 


DU ші 
2 
for the longitudinal stretching. A trignometric identity 
is employed yielding 
X = x + A(1 - cos kx) 
This research added the capability to control the location 
of the refinement through 
X = x + A(l - cos(kx 4 6)) 


The map factor 23 





15 defined by 


ax = | 
Е: = 1 + kA sin(kx + 6) 


78 





The maximum and minimum values are 


ЭХ 14M, X = 1 - Ka 
max “min 


The ratio R of maximum map factor to minimum map factor is 





_ I+kA 
re 
Solving for A yields 
2 В-1 
A= X(RH 


For purposes of this experiment R was chosen to vary from 
1.0 to 4.0. A value of R= 4 implies that the minimum X is 
one fourth of the maximum X. The Y transformation 15 
performed ina similar fashion except that Y = y * B sin Ly 
is employed where 


Y = the transformed grid, 


y = the original grid, 
B = a constant, and 
L = 2n/W 


Placement of the high resolution is accomplished by 
determining the value of 6 required to place the minimum map 
factor as desired. For instance, if the refined grid is 
desired in the middle of the channel (L/2) then 5 would be 


equal to 1/2. 
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Figure 22 is the 48-h forecast with the triangular 
Subdivision for a uniform grid. This forecast will act as 
the control for the triangular cases. Figure 23 is the 
triangular subdivision, 48-h forecast with R = 2.0. Figure 
24 is the triangular subdivision, 48-h forecast with 
R = 3.0. Figure 25 is the triangular subdivision, 48-h 
forecast with R = 4.0. All of the 48-h u and v fields 
appear identical. This indicates that there is no degra- 
dation in forecast skill when using the same number of 
degrees of freedom and locating many of those unknowns into 
a refined area. There is no noise apparent in the 
transition regime. However, there are definite differences 
in the divergence field, although the magnitude of the 
maximum/ minimum value of divergence is relatively small 
055107 7 

Table 4 shows the harmonic analysis amplitudes of the 
geopotential fields for the different triangular cases. The 
relationship between initial and final amplitudes remains 
the same regardless of the degree of variability. The phase 
propagation of the v field in the control case is 97.7%, 
whereas it is 96./%, 95.4% and 94.1% for R values of 2, 3 
and 4. The variation from the control of the phase propaga- 


tion with grid refinement is under four percent. 
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Fig. 24. As in Fig. 22 for R = 3.0. 


83 


0 





Y-AXIS (METERS( »10° 


Y-AXIS (METERS) «10" 


x10" 


7.5 


Y-RXIS (METERSI 


7.6 


6.5 


5.5 


4.5 


7.5 


5.5 


4.5 


6.5 


5.5 


4.5 


(СВЕ ЕМЕ ТОАТ 
TAU - 49 


12122115 


Z 


0.0 


Y 


1.0 






FS 


2.0 3.0 4.0 .0 
X-AXIS (METERS) мо 





ща 


5 


их 


VORTICITY 
TAU = 48 


2.0 3.0 4.0 
X-AXIS (METERS) 


29. 


12122115 





S.0 
x10 


As in Fig. 


84 


22 


SE Dd NES 
TAU = 48 





12122115 
a 
om 
>< 
=< 
x 
e 
1.0 13.0 
X-AXIS 
TAU = 48 
12122115 
bal x 
тэ 
9 
5/8 | 
эж 
š [n 
ael | 
x" 
T 
> | 
5; 
0.0 1.0 2.0 3.0 4.0 5.0 
X-AXIS (METERS) міс“ 
TAU = 48 
12122115 


«10 
7.5 


6.5 


Y-AXIS (METERS) 
5.5 





1.0 


2.0 3.0 4.0 
X-RXIS (METERS) 


for R ding. 





TABLE 4 
Harmonic Analysis (Geopotential m/s?) 
Triangles 
R = 1.0 R = 2.0 R = 3.0 R = 4.0 
Wave Initial 48-H Initial 48-H Initial 48-H Initial 48-H 
1 гт 64.65 67.91 64.95 67.91 63.70 
„00 455 “00 „69 „00 „99 200 57 
нэг 72.55 13.41 13.24*13.29 11.16 


> Cc no 


.00 „21 „00 2 .00 „54 . 00 . 28 
5 1.90 1.70 №267 1930 1.39 1-6 13 2 22. 
6 „00 .16 .00 ЁЗ .00 136 „00 22] 


Figure 26 is the control case for the rectangular 
version, 48-h forecast. Figure 27 is the rectangular 
subdivision, 48-h forecast with R = 2.0. Figure 28 is the 
rectangular subdivision, 48-h forecast with R = 3.0. 

Figure 29 is the rectangular subdivision, 48-h forecast 

with R = 4.0. Close inspection reveals no noise in either 
the % u or у fields. No degradation in forecast skill has 
been observed due to increased resolution with a rectangular 
Subdivision. All Ф, u and v fields are very similar in 
structure. The main difference in the forecasts lies in the 
divergence field where the magnitudes are small ( 1078571) 


because of the specified mean depth. 
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Table 5 is similar to Table 4 except for the 
rectangular cases. Phase propagation for the control case 
у field is 98.0%, whereas it is 98.3%, 98.1%, and 9/.9% for 
R values of 2, 3, and 4. The variation from the control 
case of the phase propagation with grid refinement is less 
than half a percent. 
TABLE 5 
Harmonic Analysis (Geopotential m/s?) 
Rectangles 
R = 1.0 R = 2.0 R = 3.0 R = 4.0 
Wave Initial 48-H Initial 48-H Initial 48-H Initial 48-H 
1 68.03 67.48 67.98 65.26 68.22 66.45 68.14 66.34 
20.1 1.00 0:3 . 36 .00 285 -00 122 
ИГ 12.05 13.39 15726 713.49 13.81 13.36 12.39 


> c N 


.01 „46 „05 вой Koi „48 „00 „43 
5 ШЕ): a8» lw5Z 09 4.54 -1:59 1.46 1.19 
6 .02 25 „07 „04 .01 .30 .00 #22 


The tests thus far indicate a successful grid refine- 
ment capability for both triangles and rectangles. There 
are neither drastic impairments nor improvements in either 
phase propagation or amplitude changes for either technique. 
Neither version stands out above the other although the 
rectangular version has less phase propagation variation and 
better amplitude conservation. Both versions allow compara- 


ble grid refinement through the use of the procedure 
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previously described. Both allow an easy and straightfor- 
ward implementation of the procedure. 

However, there is a more fundamental difference between 
the versions. Operators, such as a simpie derivative, have 
weighting coefficients associated with finite element appli- 
cations as well as with finite difference schemes. For 
instance, a simple-centered, one-dimensional fin' te differ- 
ence derivative has weights of 1/2 and -1/2. The finite 
element weights are identical to the finite difference 
weights in this one-dimensional case. The weights can 
similarly be computed for two-dimensional operators such as 
a Laplacian. The quadrature formulas used to determine the 
weights are explained in more depth in Apperdix D. The 
geometry of the triangles affects the weights. For 
instance, the two-dimensional Laplacian applied to a finite 
element application on an equilateral triangular subdivision 
gives weights of 1/6 to all surrounding points, and -1 to 
the center point, as indicated in Figure 30. However, the 
weights change if a triangle with aheight equal to one half 
the base is employed as shown in Figure 31. The fact that 
some weights are zero is a strong cause for alarm. If the 
triangles are flattened further, the farthest end point 
weights become negative. This is strictly a function of the 
geometry of the triangle and isa clear sign that use of 


the triangles warrants extreme caution. In fact, Williams 
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Fig. 30. Weights associated with a Laplacian operator 
for a triangular subdivision using equilateral 
triangles. 
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Fig. 31. Weights associated with a Laplacian operator for a 
triangular subdivision using triangles with a base 
equal to twice the height. 





(unpublished notes) has determined that the geometry of the 
triangles dictates that unless a triangle with at least two 
equal sides (isosceles) is employed the boundary conditions 


will not be exactly satisfied. 


Б. EXPERIMENT 3 

The capability of the GFEM model to resolve atmospheric 
phenomena near the scale of the smallest grid length will be 
demonstrated. This experiment is in response to Task II of 
the hypothesis. 

The ultimate motivation for any “improvement” in a 
model is a better forecast. During tne past decade, 
improvements have come, to a certain degree, by increased 
resolution. While it is intuitively obvious that increased 
resolution will improve a forecast, demonstrations nave not 
been made to show how the model resolves features near the 
smallest grid length. Staniforth and Mitchel! (1978) 
increased the resolution with a variable grid out resolved a 
synoptic scale feature in a 37x37 uniform fine mesh area. 

In this experiment, a source term simulating a mass 
source is added to the continuity equation. An expression 
for the source term in terms of a wave form is derived by 
theoretical means. The wave form allows avery small scale 
feature to be simulated. A detailed description of this 


procedure can be found in Chapter IV-D-2. By assuming the 
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wave is in equilibrium and steady state, the solution and 
source are known for all time when the Rossby number is 
small. 

The development of the source term expression followed 
traditional quasi-Jeostrophic theory. One requirement for 
quasi-geostrophic theory is that the Rossby number be small. 
In the linear case, the Rossby number will be approximately 
0.1. However, in the nonlinear case the Rossby number is 
approximately 0.4 and the quasi-geostrophic theory 
assumption is violated. 

The first test will be on a uniform grid for the 
rectangular GFEM and finite difference models and the second 
test will be on a variable rectangular grid GFEM model. 

High resolution versions of the uniform test will be run and 
act as the control. A mean depth of 5 km was selected. 
Perturbations of 2.5 m/s and 25 m/s upon a mean flow of 

10 m/s wili be evaluated for the uniform grid. This will 
allow both the linear and nonlinear aspects of the different 
models to be observed. Here linearity is implied by the 
smallness of the perturbation amplitude. The source term 
solution is nonlinear and nonlinear interactions do occur 
when the amplitude is small. Each test will be integrated 
for 96 h of forecast time. A perturbation of 25 m/s upon a 
mean flow of 10 m/s will also be evaluated for the variable 


grid. 
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The main thrust of this experiment is to determine the 
grid-scale sensitivity for the GFEM model and the equivalent 
finite difference model. Sensitivity will be measured by 
insertion of a source term into the continuity equation. 
However, the theoretical development presented in Chapter 
IV-D-2 is for the rotational component of the wind. Because 
the initial conditions for u and v are non-divergent, the 
divergence in the model must increase from a zero initial 
state. This increase will require an adjustment process. 
This serendipity effect will be exploited during the 
comparisons. 

The linear cases examined represent a small-scale short 
wave embedded in a long wave pattern. Using a small per- 
turbation for the first tests will establish a level of 
confidence that each model responds to the source term 
creditably. The nonlinear cases represent an active vortex 
in amean flow. The maximum u component is 35 m/s whichis 
hurricane/typhoon strength. Physically, the source term in 
the divergence field opposes the advection of the vorticity 
field, and “anchors” the steady-state vorticity solution. 

In either the linear or nonlinear tests, the source 
term is nonlinear. By constraining the perturbation to be 
snall, the source term plays a linear role. However, when 
the perturbation is large, the source term allows full 


nonlinear interactions. 
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Е. RESULTS 3 
ша Linear Case - Uniform Grids 

Figures 32 and 33 are the initial and 96-h fore- 
casts for the rectangular GFEM model. There are 156 degrees 
of freedom in these forecasts. The initial v field has a 
maximum value of *2.27 m/s across a separation of 4 incre- 
ments. The initial u field maximum and minimum are 12.5 m/s 
and 7.5 m/s. The source term is shown in the initial fields 
and remains constant throughout the integration. The diver- 
gence value increased to a steady state solution during the 
first 12 hours with an oscillation that died out after hour 
24. Figure 34 is the equivalent finite difference 96-h 
forecast. Figure 34 demonstrates that the finite difference 
model performs well for these initial conditions. Figure 35 
15 ane high-resolution 96-h forecast for the rectangular 
GFEM model. Close inspecticn of Figures 33 and 35 shows 
little difference. The low resolution forecast with 156 
degrees of freedom has converged to the high resolution 
forecast with 600 degrees of freedom. Figure 36 is a graph 
of v component amplitudes at hour 96 as a function of x and 
y wave numbers for the low and high resolution and finite 
difference models. The high resolution (8253115) and low 
resolution (S153115) forecasts are in excellent agreement. 


However, the FDM (F153115) forecast departs from the control 
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Fig. 33. As in Fig. 32 for a 96-h forecast. 
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even with this relatively linear case. For uniformity of 
comparison, similar ranges of x, y wave numbers will be 
shown in Figures 36 and 41. 

2: Nonlinear Case - Variable Grids 

Figures 37 and 38 are the initial and 96-h 

forecasts for the rectangular GFEM model with low 
resolution. The initial v field has a maximum of *22.7 m/s 
across a separation of 4 increments, but most of the change 
occurs over 2 increments. The initial u field maximum/ 
minimum is 35 m/s and -15 m/s. The source term as shown 
remains constant during the integrations. The adjustment 
process is evident when viewirg 96-h forecasts. The final v 
field has a maximum of 445.3 m/s, and a minimum of -38.0 m/s. 
The final u field has a maximum of 40.2 m/s and a minimum of 
-29.8 m/s. Figure 39 is tne equivalent finite difference 
model forecast of Figure 38 and the marked differences in 
the forecast are readily apparent. The maximum final v 
component is 41.2 m/s and the minimum is -47.3 m/s. The 
maximum final u component is 50.6 m/s and the minimum value 
is -38.6 m/s. Figure 4G is the 96-h forecast for the 
rectangular high-resolution GFEM model. Comparisons of 
Figures 38, 39 and 40 show a convergence of the low- 
resolution GFEM version towards the high-resolution solution 
with a much poorer showing for the finite difference model. 
Figure 41 is a graph of v component amplitudes at hour 96 as 


a function of x and y wave numbers for the low and high 
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V component amplitude as a function of x and y 
wave number at hour 96 for the high (5258115) and 
low (R158115) resolution GFEM models and the 
finite difference (F158115) model. Perturba- 
tion = 25.0 m/s. Contour interval is 0.5 m/s. 





resolution and finite difference model. As in the linear 
case, there is strong agreement between the low and high 
resolution forecasts. The departure of the FDM from the 
control is readily apparent in this highly active case. The 
inability of the FDM to properly handle the nonlinear 
interactions manifest itself as large spurious amplitudes at 
high x and y wave numbers. 

The relatively poor showing of the finite difference 
model for the forced cases shown so far can be improved, to 
a certain degree, by increasing the resolution. However, 
increased resolution requires more computational effort. 
Figure 42 is the 96-h forecast with a perturbation of 25.0 
m/s for the FDM model. The initial conditions are identical 
tc the nonlinear cases shown in Figure 37. The difference 
in the forecasts is the resolution where there are now 576 
degrees of freedom (24X24). The 96-h forecast is certainly 
better than the 12X12 forecast. However, there is some high 
frequency noise in the forecast. Figure 43 is identical to 
Figure 42 except that there are 1296 degrees of freedom 
(36X36). The high frequency noise is readily apparent and 
1s a manifestation of nonlinear aliasing. For this forced 
case the finite difference model is unable to achieve the 
Same forecast as the low resolution GFEM model. This result 
will be amplified when a model with the same number of 
degrees of freedom as the low resolution model is employed 


but with a variable grid. 


108 


Y" une u e 1 


= 
1 
а 
7 
1 





ED 


a 
E 
froin 
Б” > 
= u, 
E 


FS25 


SPRER 


GEOPOTENTIAL HEIGHT 
05258118 


0792 


КС 





25.0 


C 
Ze 


E 


X-RXIS 
Mus 35 


— 


1.0 


БІХН-А 


"10° 


Sha: 


4.0 


3.0 
“AXIS (METERS) 


WN 


АА 


ё 
К) 


ХҮС 22772 
N 2 


< Z 
SIEH 


2 


4.0 


X 


1.0 


0.0 


ЕСУ 6% S£ 


БЕ 5 
(SYI13H) S]XU-A 


U 


TAU = 96 


3 


325811 


5258115 


хуижххикхумим мина дижхкк 


m Aa ISP AA A ARRAK 


5 
wit 


3.0 4.0 


2.0 
Х-НХ15 (METERS) 


DIVERGENCE 


1.0 


X X XX X x X XXX X X X X X 3 
X X x X X X X X X.X X x< 


X X X xX X X Xx 
X x X X X X X X W X 





0.0 


са ss 5% S'S 6% бв 
Kir (SM3JI3H) SIXU-X 





XXXXXXXXXXEXXXXXXXXXXX 
X XX X X X X X X XXX X X X X X X X X X X x 
Xr ON oos хххххх 


5.9 






x x KXXHX XxX XXXXX 
Xx XXX 









4.0 


X-RXIS (METERS) 


үш И ЛҮ 





Az SH = 
25 м, = 


Жос 
у f SUN ^ EY x 
ЧА (ТН ) хўх 
220822 : 
S Er 









3.0 


ае ны 


ИН ohn MX x x 
4.11 k yx 
ONE X X X x Oe x< 
i x X XXX x 


XF x X xx xxx ххх 
YX XX AMEX XX x XXX 
хх хх X X x X X x 
X X X X ⁄X xxx X X X X X X X X X X X 






2.0 






0 






1. 






g 





s'g sé 8% S £ 


579 5:5 
(643134) S1XU-A 


“10 


20115 


ТШ = ЧЕ 
FSZ 


TRU + "96 
3258115 





ххх хххх > 
АВ A N A X X K X X оо 
' y 

o 
= m 
E 
ob 
"mx 
u 
2х 
ag 
>< 

Q 

e 

e 

578 SL 40 55 су е“ 
‚Dix ISY3LMI SIXH-A 

eb 
x 

o 
- 
у 
о: 
ы 
sb 
nz 
у 
x 
IE 
r i 
>< 

X x x x x 

X X X X X x X X o 

- 

о 

S 





85 5: 65% 
(CH3134) CIXU-A 


39 with 576 degrees of freedom. 
109 


As in Fig. 


42. 


ке 





Figs 


GEOPOTENTIRL HEIGHT 
ТАУ = 96 


5356115 


OO OCC ie) © 







uw Оба 
N 122245 AR IRES SHOT TRIES $9328 3 9 90 obs 
> 20509750 KA BO SSS SF EGE 
да 33 SEE EET 
2 
в НН 
Ош | го 
ә 
НОЯ 
ен I TIE 
о EE Pee 5955555 
"BUS PNIS 
Е НАЕ d 8252255150 525 : 
” оо Ak 59% 
be] е A 
0.0 1.0 $0 —— 3.0 4.0 S.0 
X-RXIS (METERS) мб 


ТРО - 96 
FS358115 


ма” 


8.5 















EUN RARE SE ETE IEEE 
я 252255221010) 252555 


9 $06 X 25 


22 
4994 ? 2 


7.5 


уа 
RESTE 


6.5 


HE AAA 
999990 


30202099 
ere 


5.6 


ОВО 
2СХУуО0О009о00099О0О099 
9 229099 


Y-AXIS (METERS) 


4.5 


оз; Ф е 
ООС м осо Оо ООО 
(445252076 ets Y 


3.5 





I 1.0 2.0 3.0 = e 5.0 
X-AXIS (METERS) “С 
ШЕШЕСІ I Dy. 
ТАУ - 96 
е 5158115 
llo 
хе! 

к ë 
и ка 
Би! Shu 
74 а: 
= e 
Qu! 4 
27 MESH 
» | Хээ 

s a 

w EX 426 

^ ER | - 

0.0 1.0 2.0 0 


3.0 4. 5.0 
X-RXIS (METERS) "10° 


Y-AXIS 


мо! 


Y-AXIS (METERS) 


x1G 


Y-AXIS (METERS) 





STREAMLINES 
TRU = 96 


FS35811S 





X-RXIS 


TAU = 96 
5358115 





2909 ОО ty 
o 2099999944 206009442039400о09О00090090 


Ф „ 
5258294 0 006056 
” ` 

Е 


200909 


Ф 
99220оЭ922о9292429454096 МО 


679 020212 
2) 


Во 10 -20 3.0 0 5.0 
X-AXIS (METERS) x10' 


DIVERGENCE 
18358115 


aS pe ee yt a кх КАҢ, 
EEE See ce IA 
25 ee 29 22355552 x 
A t M 


TT | 
Q S a ASE 
«Lo 4 x^ 74 
525552551725 E 
: Q а 00000 $ 4 D 


% x 9 3 50012 қ X X X 
995 Ба ы 


3.0 4.0 5.0, 
X-AXIS (METERS) “10 


> 


43. As in Fig. 39 with 1296 degrees of freedom. 





This comparison of the finite element and finite 
difference methods for the source term dramatically high- 
lights the potential improvement when utilizing a Galerkin 
approach. Both graphicai displays and harmonic analysis 
leave little doubt as to the superiority of the GFEM model 
over the equivalent finite ditference model. 

Adding an extra margin to the already established 
superiority can be accomplished by using the low-resolution 
model with a variable grid. Figure 44 is similar to the 
tests in Figure 38 except that a variable resolution grid 
(В = 2.0) was utilized. A 96-h forecast for a variable 
resolution grid (R = 3.0) is shown in Figure 45. Comparison 
of Figures 38, 40, 44, and 45 show tie convergence of the 
low-resolution uniform solution towards the high resolution 
solution as the degree of variability is increased. Figure 
46 isa graph of geopotential amplitudes at hour 96 versus y 
mode wave number (x wave number one) for the uniform, 
variable, and high resolution models. The R = 2.0 case 
better represents the control than the uniform low 
resolution forecast especially at low wave numbers. 

However, at wave number three and above all models have 
Similar differences from the control. 

The harmonic analysis demonstrates an improvement with 
increasing variability but is not conclusive. Another 
method to evaluate the improvement with increased varia- 


bility is to look at the difference charts between the 
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control and each level of variability. The difference chart 
in Figure 47 is the contro! 96-h forecast (high resolution) 
minus the low-resoiution 95-h forecast. Figure 47 shows a 
maximum difference ІП the center of the vortex. Figure 48 
is the comparison of the high resolution 96-h forecast and 
the R = 2.0 variable grid S6-h forecast. The maximum 
difference has decreased from the uniform case. Figure 49 
is the comparison of the high resolution 96-h forecast and 
the R = 3.0 variable grid 96-h forecast. The difference 
centers at the right-hand edge in Figures 4/7, 48 and 49 are 
a result of post processing interpolation errors. Although 
the maximum difference in the center of the vortex is 
greater in Figure 49 than figure 48, the latitudinal 
difference gradient is less. Root mean square differences 
from the high resolution control are 652.2 п2/52, 609.8 
m/s? and 469.7 m/s? for the low resolution, R = 2.0 and 

R = 3.0 tests respective:y. The forecasts have definitely 


been improved at each level of increased resolution. 


G. EXPERIMENT 4 

This experiment wil! demonstrate a moving grid 
capability. Harrison (1973) demonstrated the use of a 
moving grid in his finite difference model. The inherent 
problems with moving finite difference models are abrupt 
Changes in resolution causing spurious noise. The noise, 


unless damped, is magnified if left unattended every time a 


ІШЕ» 





ка. 


47. 


GEOPOTENTIRL HEIGHT 
ІН) = 96 


5258115 





GEOPOTENTIRL HEIGHT 
TRU = 96 


5159115 


Y-AXIS 





DIFFERENCE CHART 
IRU = 96 


258 - 158 





Forecasts and difference charts at hour 96 for the 
uniform high (5258115) and low (151598115) résolu- 
tion GFEM models. Perturbation = 25.0 m/s. 
Contour interval is 600 m2/s® for the geopotential 
fields and 200 m ЁС for the difference field. 


116 





БЕШ ШШЕ HEIGHT 
TAU = 36 


4 2 
5258115 





GEOPOTENT TAL HEIGHT 


S178115 


Y-AXIS 





DIFFERENCE CHART 
ТНУ « 96 


258 - 178 





Fig. 48. As in Fig. 47 for the hi | 
17. gh resolution (S258115 
and the variable (S178115) (R = 2.0) udi ца... 


117 








GEOPOTENTIAL HEIGHT 
TAU = 96 


5258115 


21.0 


16.0 


Y-AXIS 


11.0 





1.0 6.0 Е 16.0 21.0 
GEOPOTENTIRL HEIGHT 
IAU =- 96 
$198115 


Y-RXIS 
16.0 21.0 


11.0 


6.0 





16.0 21.0 


сте 
DESEE SIE ORT 
288, = 138 





Fig. 49. As in Fig. 47 for the high resolution (S258115) 
and the variable (S198115) (R = 3.0) GFEM models. 


118 





grid moves. In a GFEM model there is a smooth transition 
into the fine grid area with no attendant noise. Therefore, 
moving the entire grid with the fine grid is simply a matter 
of accurate interpolation. 

Three tests will be performed. First, uniform low and 
high resolution models will be integrated. These models 
will not use a moving grid. Then amodel with the same 
degrees of freedom as the low-resolution model, but 
employing a moving, variable grid, will be run. The high- 
resolution model will act as the control. Identical initial 
conditions will be used for ali three tests. These condi- 
tions are similar to those developed for experiment 3, 
except that the source term wi!l be zero throughout the 
course of the integration. The 1п1:1а1 conditions include a 
sharp trough with a relative vorticity maximum which is 
advected by the mean flow in the absence of any forcing (the 
source). 

Of particular concern is the generation of noise 
incurred during grid movement. Therefore, a forecast of 
96-h is made during which this grid moved eleven times. As 
previously stated, the model utilizes absolutely no diffu- 
sion, filters or friction. If spurious modes are excited 
during the movement of the grid, they will appear in the 


harmonic analysis. 
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The procedure of moving the finite element grid is 
achieved by displacing all the nodes one unit distance in 
the direction desired. This retains the relative 
distribution of the nodes and merely moves the entire grid 
as an entity. Ramifications of this approach are: 

live The areas of each element remain the same and the 

numerical quadrature coefficients for all the 


associated integrals do not have to be recomputed. 


2 The unit distance moves should be the smallest 
grid increment; and 


Ве: Direct solvers that utilize a preprocessing 
procedure based on the grid distribution will not 
require the preprocessing to be recomputed after 
each move. 

The most important aspect of moving the grid is to accurate- 
ly interpolate the new field values for the new grid point 
position. In this dissertation a standard International 
Mathematics and Statistical Library's (IMSL) bi-cubic spline 
interpolating routine was employed. The movement will occur 
only in the x direction and, therefore, a one-dimensional 
interpolator suffices. The interpolator accommodates the 
periodic boundary condition. However, interpolators for 
non-periodic boundary conditions are available. The accu- 
racy of the interpolator is exact at grid points for a 
uniform distribution. Initial experimentation verified that 
a field could be moved with no loss in accuracy for a 


uniform grid. Experiment 4 will test its capability to move 


a variable grid with R = 2.0. 
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Н. RESULTS 4 

Figures 50 and 51 are the initial and 96-h forecasts 
for the uniform low-resolution model. The perturbation is 
25.0 m/s, mean flow 10 m/s and mean depth 5 km. The initial 
divergence is a zero field but spins up through geostrophic 
adjustment. The northeast/southwest trough orientation is a 
manifestation of the nonlinear interactions occurr:ng for 
this highly active perturbation. Figure 52 is the 96-h 
forecast for the uniform high-resolution model. The 
northeast/southwest trough orientation is evident in this 
forecast as well. Figure 53 is the 96-h forecast “or the 
low-resolution, moving variable grid (R = 2.0) model. 

Time steps were 4320s, 3085.7s and 3600s for :he low, 
high and variable models, respectively. The variable model 
took 95 full time steps and two Mir steps during the 96-h 
forecast. The half steps are required on!y to start the 
model to obtain the N and №+1 time levels in the leapfrog 
time stepping. The grid moved 11 times during the forecast 
at time steps 8, 16, 24, 32, 40, 48, 56, 64, 72, 80 and 88. 

The x axis origin in Figure 53 is 3530 km indicating 
eastward grid movement. The vortex is still in the fine 
mesh area. Additionally, there is no noise in the forecast 
fields. To reiterate a previous statement, there are no 
filters, diffusion or friction terms employed during this 


integration. Comparisons among Figures 51, 52 and 53 show 
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that tne variable grid forecast more closely approaches the 
control than does the uniform low-resolution model. This is 
to be expected since the fine mesh area has been moved to 
coincide with the active region of the forecast domain. 
Table 6 shows the geopotential harmonic analysis for 


the three models efter a 96-h forecast. 


TABLE 6 
Harmonic Analysis) (Geopotential) for Experiment 4 
Units are m¢/s¢ 


Wave Low Resolution High Resolution Moving Grid 


1 676.1 904.7 850.0 
2 167.4 152.9 146.6 
5 122.5 175-9 1159.6 
1 54.7 3910 20:10 
5 29:-3 10:22 9.6 
6 2! 2 es 


Comparisons of the higher wave number amplitudes show 
almost identicai amplitudes in y modes four, five and six 
for the moving case and high-resolution uniform case. The 
amplitude in y mode six indicates that the movement of the 
grid has not generated unwanted noise. The overall improve- 
ment for the moving case is a result of the finer grid's 
ability to better resolve the small-scale interactions 


occurring. 
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VI. SUMMARY AND RECOMMENDATIONS 


Numerical weather prediction requires that the 
equations governing atmospheric flow be solved numericaily. 
These flows cover a large range of scales. Classical 
modeling schemes have successfully forecast the longer 
scales and improvements are being sought for the smaller 
Scales. 

In this dissertation, the GFEM has been employed for 
the shallow-water equations in a channel domain. The 
differentiated form of the equations on an unstaggered 
finite element grid with a semi-implicit time scheme has 
been utilized. An equivalent finite difference model 
representing the conventional approach has also been 
utilized. Analytic initial conditions representing simpie 
atmospheric waves and a more complex source term have been 
used to test the different models. 

The hypothesis of this dissertation is that the GFEM is 
a viable option for numerical weather prediction. The 
theoreticai foundation of the GFEM leads to a minimization 
of the error between the actual equations and their approxi- 
mation. This minimization implies a basis for expected 
improvement in forecast capabilities. Three features were 
addressed to support that hypothesis. First, alternatives 


for variable grids were investigated. Second, forcing near 
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the smallest grid scale determined if the GFEM model 
resolved the smaller scales better than a finite difference 
counterpart. Finally, a proof of concept for the ability to 
move the grid was demonstrated. Each task adds to the 
usability of the GFEM and illuminates the well-rounded 
potential that it contains. 

rirst, the ease with which the method allows variable 
grids was demonstrated. In all cases, no numerical 
techniques were employed to damp or filter any field. The 
grid refinement obtained when using a rectangular 
subdivision was more attractive than with a triangular 
subdivision. It allowed refinement to be easily obtained 
with increased accuracy from the higher order polynomial and 
decreased computational effort due to reduced operation 
counts. Additionally, it was shown that the geometry of the 
triangular elements greatly influences the problem. 
Furthermore unless specific constraints are enforced, the 
boundary conditions would not be satisfied. 

Second, the ability of the GFEM to resolve atmospheric 
phenomena occurring near the scale of the smallest grid- 
length was shown. A small-scale steady state solution was 
cast into a source required for that solution. Integrations 
Showed that the GFEM model resolved the small-scale forcing 
admirably and far exceeded the performance of an equivalent 


finite difference model. Even when the resolution of the 
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finite difference model was doubled and tripled, it could 
not match the performance of the low-resolution variable 
grid GFEM model. The responsiveness of the GFEM model near 
the smallest resolvable grid length clearly demonstrates 
that the variable grid concept is worth the computational 
effort. 

Finally, the concept of moving the variable grid was 
demonstrated. The ability to resolve a small-scale 
phenomenon is a tremendous asset to the GFEM. The ability 
to move the grid with the atmospheric phenomena further 
enhances its usefuiness. 

Througnout this dissertation the GFEM model has 
demonstrated desirable characteristics. It has conservative 
properties. It propagates atmospneric waves better than an 
equivalent finite difference model. It allows variable- 
resolution grids and responds better than an equivalent 
finite difference model near the smallest gridlength. 
Moving grids can be achieved with no apparent noise 
generation. It can utilize direct solvers and is a natural 
choice for vectorization on iarge computers. It truly is a 
viable option for simulation of atmospheric flows. 

The success of this dissertation mandates further 
research. The next logical step would be a baroclinic 
model. Completion of a baroclinic model should allow its 


application as a regional model. Additionally, different 
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basis functions could be employed for the vertical 
discretization, such as solutions to the vertical structure 
equation. For the moving variable-resolution grid, further 
research into two-dimensional movement is required. The 
proof of concept in this dissertation of one-dimensional 
movement can logically be extended to two dimensions. Upon 
completion of two-dimensional movement, an Operational 
forecast capability for tropical cyclones should be investi- 
gated. The superior small-scale response of the GFEM 
indicates potential increase in skill for both regional and 


tropical cyclone/hurricane prediction. 
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APPENDIX A 
STABILITY ANALYSIS 


A one-dimensional linear analysis of the forecast 
equations will be presented. The analysis was obtained 
through personal communication with Dr. Andrew Staniforth 
(Recherche en Prevision Numerique). The primitive form of 
the forecast equations and a semi-implicit time scheme wil! 
be analyzed, because the results are identical for the 
vorticity/divergence form. The one-dimensional equations 


with a mean flow U are 


ди ав edi + ` 
3t эх лү vor 
ду 27 рим _ x 
3t Ec fu (A-2) 
9% 94 22129 

5t эх Ú sx (A-3) 


Time derivatives are evaluated with a centered time dif- 
ferencing and the other terms on the left-hand side are 
averaged between time levels (t + At) and (t - At). 


Equations (A-1), (A-2) and (A-3) become 


u(x,ttAt) - u(x,t-at) E l Me lxtax, ttat) - d(X-Ax, ttAt) 
2 


2At 2AX 


САХ 


+ (9(xtAx, t-At) - $ (x-Ax,t-At) " m rU(XtAX,t) - u(x-4Ax,t 
15--11 PAX ] 


шиг) (А-4) 
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v(x,ttAt) - v(x,t-4t) . jj [vorax t - у(х-Ах,% 1 
2At 2АХ 


- f u(x,t) 


o(x,ttat) - o(x,t-at) , 


eat 


гою 


[Au GxrAx, trat) - u(x-Ax,t*At) 
2AX 


» (u(xtax,t-At) - u(x-Ax,t-At 3 .. D($(<tax,t) - $(x-ax,t); 


2AX 2AX 


Assuming that a function, F, can be expressed as 


Ext) = p'el(kx +wt) 


then 
” i(kx + wt) 
ҰНЫ = Ex taat) Е = sin(wAt)e 
l i(kx + wt) 
F(x,t+at) : F(x,t-^t) . F' cos(wat)e 


Equations (A-4), (A-5), and (A-6) become, after 


substitution, 


isu - fy! + ikeg' = 0 
l isv' = 

kurt ТЕ 0 

T ! 159” = 

ikcou Ет 0 


(A-5) 


(12) 


(А-7) 


(А-10) 


(A-11) 


(A-12) 





Where sessinwwAt * k'uAt, c = cos wAt and k' = sin ——, 
The determinant of the system of equations (A-10), (A-11) 


and (A-12) when set equal to zero yields 


2 
S$ rs 2 22 B (A-13) 
дъ Со (Е + Ксе)1 = 0 
The roots are s=0 and 
s? - (fat)^ * c^(kat)? à (A-14) 


The roots of this second order equation when requiring w to 


be real yield the stability criterion 


Т. ааа (А-15) 
1-1 6 
AX 
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APPENDIX B 
NUMERICAL QUADRATURE 


Evaluation of the Galerkin integrals requires a fast, 
efficient method . Two different elements, triangles and 
rectangles, have been employed and the quadrature schemes 
for evaluation are presented nere. In the first case, 
triangles, area coordinates are used. For the rectangles, 
integration formulas are based on an orthogonal axis trans- 
formation. In both cases, the integrals to be evaluated 
contain either products of basis functions, products of 
derivatives of basis functions, or a mixture of both. 

Zienkiewicz (1971) daveloped the relationship between 
the Cartesian coordinates of a triangle and the area 


coordinates as 


(B-1) 
KL, + LK, + Lak, 
Y = L,Y, + Lavy * 374 (8-2) 
1 = Li + L; + L; 

(B-3) 


where (x1,41)» (X2,y2) and (x3,y3) are the Cartesian coordinates 
сие пещи още 5” мегстсез ад Lj, Lo and L3 are the area 


coordinates. Here 


Е A,/A 3 = А /А 3 = Az/A 
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where A), A, and À4 are shown in Figure 54 and A is the area of 


the entire triangle. He further shows that 


f[ 4 1221,5 dxdy = tay 2А 
The L's are the basis and test functions used during the 
approximating process. Figure 55 defines some necessary 
parameters for evaluation of the derivatives. 


Equations B-1, B-2 and B-3 can be written in matrix 


form as 

[A 2A b. а, Ш 

m 
La 2A р, а; Y 
Differentiation of Eq. (B-5) shows that 
3 bD. 

9 1 9 
L<= § — —- (B-6) 
9 : 3; 9 (B {| 
Е > — — = 
9Ё 1 2A oL. 


However, the derivative of the basis function with respect 
to the area coordinate is non-zero only where the basis 
function and area coordinate coincide. In fact, the non- 


zero value is one. Therefore, using Y as a basis function 
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Fig. 54. Cartesian coordinates vs. natural coordinates. 





Fig. 55. Triangle definitions for area coordinates. 
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All integrals can be evaluated by combination of these 
quadrature formulas for triangles. 

An orthogonal axis transformation will allow similar 
quadrature formulas for rectangles. The rectangle shown in 


Figure 56 can be transformed by 








aS a nm b (B-9) 


The values of с and n at each corner are shown in 


parentheses. The basis function, V;, becomes 


yo, tege m) (8-10) 
i 4 





Fig. 56. Orthogonal axis transformation for rectangular 


integration formulas. 


Derivatives of the basis function become 


E UE 
ЭХ а эс 
апа 
Aei (6511) 
ду b әп 
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By substituting and integrating, the 4 x 4 matrix with the 
interaction coefficients can be determined for a derivative 
or straight inner product. For instance, the straight inner 


product is 


TAI 
МЕ =] ч xis = аЬ / | шэг dzdn 


1 


= 


II 
— 
e| 
— 
го 
со, 
уч 
JN 
ши 
on 
NO 
«deo 


(B-12) 


and the mixed derivative is 


my 
di 1 -ff x dxdy = 30 [= —À. dzdn 


1 
-f уе | (180, 11131) 41 


р 2 2 


These quadrature rules allow the evaluation of any integrals 


when using rectangles. 
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